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AN  EFFICIENCY  STUDY  OF  SEVERAL  TECHNIQUES  FOR  THE 
NUMERICAL  INTEGRATION  OF  THE  EQUATIONS  OF 
MOTION  FOR  MISSILES  AND  SHELL 

ABSTRACT 

The  numerical  integration  of  the  equations  of  motion  for  missiles 
and  shell  is  a  frequently  occurring  computational  problem  in  government 
laboratories  and  in  the  aerospace  industry.  Because  of  the  repetitive 
nature  of  such  computations  and  the ■  continuing  requirement,  it  is  desir¬ 
able  to  make  these  computations  as  efficiently  as  possible  in  order  to 
minimize  the  computer  time  spent  in  their  solution.  The  cost  in  terms 
of  time  required  to  solve  representative  problems  to  a  specified  accuracy 
is  a  measure  of  the  relative  efficiency  of  numerical  integration  methods 
and  can  be  determined  only  by  experimentation.  This  report  describes 
the  results  of  such  an  experimental  study  and  finds  that  the  Kutta-Merson 
procedure  with  automatic  and  continuous  interval  adjustment  is  far 
superior  to  the  predictor-corrector  techniques. 
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INTRODUCTION 


Prior  to  the  study  culminating  in  this  report  the  Ballistic  Research 
Laboratories  Electronic  Scientific  Computer  (BRLESC)  had  available  in 
its  compiler  two  subroutines  for  the  numerical  integration  of  a  system 
of  ordinary  differential  equations.  These  subroutines  were  RKG  and 
RKGMA.  The  RKG  routine  is  the  Runge-Kutta-Gill  procedure  as  described 
by  Romanelli7*.  RKGMA  is  a  procedure  using  the  Adams  fifth  order 
predictor  formula  with  no  subsequent  correction.  The  procedure  uses  the 
RKG  routine  as  a  starter,  hence  its  name,  Runge-Kutta-Gill-Modified- 
Adams**.  The  reason  for  using  only  the  predictor  cycle  in  RKGMA  as 
opposed  to  the  classical  predictor-corrector  methods  was  that  the  RKGMA 
method  requires  only  one  derivative  evaluation  per  step  and  hence  runs 
twice  as  fast  as  the  conventional  method  when  both  methods  use  the  same 
step  size.  This  allows  a  step  size  of  h/2  for  the  same  cost  as  the  con¬ 
ventional  method  using  step  size  h  and  results  in  a  technique  with  greater 
accuracy.  This  observation  was  "verified"  in  a  study  performed  in  the 
Computing  Laboratory  and  led  to  the  incorporation  of  RKGMA  into  the 
BRLESC  compiler.  Of  course  all  this  neglects  the  very  important  consid¬ 
eration  of  stability  and  in  fact  Hamming3  shows  that  for  a  simple  predictor 
formula,  such  an  application  can  be  unstable. 

The  use  of  RKGMA  in  the  integration  cf  trajectories  nevertheless 
produced  a  very  accurate  method  and  was  initially  thought  to  be  more 


♦Superscript  numbers  denote  references  which  may  be  found  on  page  37. 
♦♦Hildebrand3  refers  to  the  Adams  fifth  order  predictor-corrector 
method  with  iteration  as  the  modified  Adam's  method.  To  avoid  possible 
confusion  this  distinction  is  noted. 


7 


economical  than  the  RKG  method.  This  conclusion  was  derived  as 
follows:  In  the  absence  of  any  practical  estimate  of  the  truncation  error 
in  the. RKG  scheme,  it  was  traditionally  used  with  a  constant  step  size 
which  from  experience  one  knew  was  adequate  and  "safe".  When  RKGMA 
became  available,  it  was  natural  to  attempt  to  solve  trajectory  problems 
with  the  same  step  size  traditionally  used  with  the  RKG  routine.  The 
result  was  that  the  method  was  stable  for  the  particular  applications 
and  yielded  a  method  almost  four  times  more  economical  than  RKG. 

Such  a  comparison  is  not  complete,  however,  in  that  one  should  also 
consider  accuracy  and  stability.  Subsequent  experiments  designed  to 
consider  these  factors  led  to  the  conclusion  that  for  trajectory  compu¬ 
tations  RKG  was  more  efficient  than  RKGMA.  These  results  were 
attributed  to  the  fact  that  to  maintain  stability  with  RKGMA  one  had  to 
use  a  smaller  step  size  which  more  than  offset  its  four  to  one  advantage 
with  regard  to  derivative  evaluations. 

A  method  for  improving  the  stability  properties  of  a  procedure 
such  as  RKGMA  arose  from  the  findings  of  Hull  and  Creemer4  who 
reported  on  a  study  of  the  efficiency  of  predictor-corrector  procedures. 
They  found  that  a  method  with  only  one  evaluation  of  the  derivative  con¬ 
sistently  seemed  the  most  efficient  type  of  predictor-corrector  technique. 
This  technique  maintained  the  advantage  of  RKGMA  with  respect  to  deriva¬ 
tive  evaluations,  but  by  modifying  the  cycling  process  improved  the 
stability  properties.  This  scheme  allows  one  to  apply  both  the  predictor 
and  corrector  for  the  price  of  one  evaluation  or  the  application  of  the 
predictor  and  then  the  corrector  step  twice  for  the  cost  of  two  evaluations. 
In  addition,  this  procedure  even  with  only  one  evaluation  provides  an 
estimate  of  the  truncation  error.  A  study  of  this  scheme,  applied  to 
the  Adams  formulas  became  the  basis  of  the  BRL  study.  An  additional 
objective  of  this  investigation  was  to  study  the  formulas  of  Crane  and 
Klopfenstein5 ,  who  derived  a  predictor  formula  which  when  used  with 
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the  Adams  fifth  order  corrector  yielded  a  method  with  increased  range 
of  stability.  During  the  study  the  Kutta-Merson  algorithm  came  to  the 
attention  of  the  author  and  was  added  to  the  study  because  of  its  favor¬ 
able  stability  properties  and  because  it  provides  a  computable  estimate 
of  the  truncation  error. 

The  studies  described  herein  were  by  no  means  exhaustive;  i.  e. , 
only  fifth  order  methods  were  examined,  but  seem  sufficiently  important 
to  warrant  reporting  since  few  such  studies  on  practical  dynamic  problems 
are  found  in  the  literature. 

The  trajectory  models  used  in  the  study  were  actual  models  currently 
being  used  by  the  Ballistic  Research  Laboratories  in  computing  Firing 
Tables  for  various  weapon  systems.  Although  the  results  of  this  study 
may  to  some  extent  be  affected  by  programming  pecularities  and  the  use 
of  a  particular  computer  it  is  felt  that  the  results  are  fairly  representative 
of  what  might  be  expected  in  general. 

General  Discussion 

Many  different  schemes  exist  for  the  numerical  solution  of 
ordinary  differential  equations.  Because  of  the  widespread  need  for 
numerical  integration,  most  computing  facilities  have  as  a  standard 
routine,  some  method  of  numerical  integration  which  can  be  utilized  as 
a  package  by  programmers.  This  leads  to  a  more  efficient  operation 
and  avoids  the  occurrence  of  redundant  programming  effort.  Because 
of  the  general  use  of  such  a  routine  or  package,  it  is  highly  desirable 
that  it  be  flexible  and  efficient.  The  flexibility  is  measured  by  the 
relative  ease  with  which  programmers  can  call  for  the  routine  and  by 
the  number  of  tasks  to  which  the  routine  can  be  applied  with  confidence. 
The  efficiency  is  best  measured  by  the  computing  time  required  to 
solve  representative  problems  to  a  given  degree  of  accuracy. 
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The  most  efficient  routines  seem  to  be  those  which  provide  an 
estimate  of  the  truncation  error  and  which  adjust  the  integration  interval 
as  the  integration  proceeds.  When  a  provision  for  interval  adjustment  is 
lacking,  the  programmer  is  forced  to  be  conservative  in  his  choice  of 
step  size  and  may  even  have  to  solve  the  same  problem  several  times 
with  different  intervals  to  assure  himself  that  the  procedure  is  stable 
for  the  particular  problem  and  the  chosen  interval. 

f 

The  truncation  error  is  a  measure  of  the  error  incurred  at  a  given 
step  and  cannot  be  used  to  determine  precisely  the  accumlated  error 
resulting  after  the  integration  has  proceeded  for  many  steps.  Despite 
this,  the  use  of  truncation  error  as  a  guide  to  choosing  the  integration 
interval  has  proven  quite  successful  in  creating  efficient  routines. 

The  computation  of  trajectories  is  a  frequently  occurring  problem 
which  requires  numerical  integration.  The  Firing  Tables  Branch  of 
the  Computing  Laboratory,  for  example,  in  1965  used  over  900  hours 
of  machine  time,  a  significant  fraction  of  which  was  used  on  the  compu¬ 
tation  of  trajectories.  These  computations  usually  allow  varying  intervals 
as  the  problem  proceeds  to  attain  a  fixed  truncation  error.  For  this 
reason  a  process  which  can  estimate  the  truncation  error  and  adjust  the 
step  size  continuously  and  automatically  can  be  much  more  efficient  than 
a  procedure  which  does  not.  The  design  of  such  a  procedure  for  general 
trajectory  computations  was  the  major  objective  of  this  study.  An  addi¬ 
tional  objective  was  the  creation  of  a  standard  subroutine  which  could  be 
used  by  programmers  in  general  on  a  variety  of  problems.  Among  the 
integration  procedures^which  have  the  most  widespread  use  are  the  two 
basic  categories: 

(1)  Runge  Kutta  methods 

(2)  Predictor-corrector  methods 
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The  general  features  which  are  discussed  extensively  in  the  literature 
pertain  mostly  to  the  following  six  areas: 

(1)  Storage  requirement 

(2)  Accuracy 

(3)  Stability 

(4)  Ease  of  starting 

( 5)  Number  of  derivative  evaluations  per  step 

(6)  Estimation  of  truncation  error 

The  Runge  Kutta  method  or  the  Gill  modification  evaluated  against 
these  features  has  become  almost  universally  adopted  as  the  best  pro¬ 
cedure  to  have  as  a  general  purpose  standard  package.  The  method  is 
very  accurate,  has  probably  the  best  stability  properties  of  the  methods 
used  extensively,  is  self  starting,  and  the  Gill  modification  of  the  method 
requires  very  little  storage  of  data.  The  method  suffers  from  the  dis¬ 
advantage  of  requiring  four  derivative  evaluations  to  advance  the  compu¬ 
tation  one  mesh  interval.  As  generally  used  the  method  also  suffers 
from  the  additional  disadvantage  of  not  providing  an  estimate  of  the  trun¬ 
cation  error  unless  used  in  a  multistep  mode1  ;  i.  e. ,  the  same  step  is 
done  twice,  once  with  an  interval  of  h  and  then  followed  by  two  steps  of 
h/2.  This  method  has  not  found  widespread  use  mainly  due  to  the 
increased  cost  of  the  multistep  application.  The  computable  estimates 
of  truncation  error  found  by  Merson6  and  more  recently  by  Scraton7 
seem  to  largely  eliminate  this  problem  however. 

The  predictor-corrector  methods  are  generally  considered  to  be 
the  most  efficient  since  in  the  conventional  application  only  two  deriva¬ 
tive  evaluations  are  required  and  because  of  the  ready  availability  of  a 
computable  estimate  of  the  truncation  error.  While  this  may  be  gener¬ 
ally  true,  there  exists  a  great  deal  of.  evidence  (this  study  for  example) 
indicating  that  for  some  classes  of  problems  the  Runge-Kutta  methods 
are  the  most  economical  if  used  with  a  technique  for  interval  adjustment. 
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Despite  the  extensive  body  of  literature  that  exists  on  the  subject , 
no  definitive  conclusions  can  be  made  when  attempting  to  choose  a  method. 
Those  studies  that  are  found  in  the  literature  usually  are  done  for  proto¬ 
type  differential  equations  and  it  is  not  generally  valid  to  extend  the  find¬ 
ings  from  such  studies  to  the  more  general  dynamic  problems  that  are 
of  greatest  interest.  Many  problems  are  such  however  that  the  most 
expeditious  solution  is  to  use  a  method  available  at  one's  computing 
facility;  i.  e. ,  even  though  the  particular  method  may  be  very  inefficient 
for  a  given  problem,  the  fact  that  it  has  already  been  programmed  and 
checked  out  offsets  other  considerations.  The  problem  may  be  a  "one-shot" 
type  which  can  be  solved  for  several  integration  intervals  without  incurring 
much  total  cost  in  computation  time.  The  trajectory  problem  is  one, 
however,  which  is  repetitive,  i.  e. ,  after  a  program  is  written  it  may  be 
used  to  compute  many  trajectories  for  a  given  projectile  or  missile  and 
quite  often  is  used  for  different  projectiles  and  missiles.  In  this  instance 
one  is  concerned  with  programs  which  are  used  almost  daily  for  a  period 
of  years  and  any  economies,  however  small,  eventually  produce  large 
accumulative  savings. 

Several  techniques  are  commonly  used  for  choosing  the  mesh  size 
in  this  type  of  problem.  One  method  is  merely  to  use  a  constant  mesh 
of  a  magnitude  which  from  experience  one  has  confidence  will  work. 

Another  technique  is  to  establish  some  empirical  rule  by  experimentation 
which  keys  on  the  magnitude  of  a  certain  variable  and  adjusts  the  step 
size  as  a  function  of  this  variable.  The  third  method  is  the  one  mentioned 
previously;  i.  e. ,  the  integration  procedure  itself  estimates  the  truncation 
error  and  adjusts  the  step  size  to  keep  the  truncation  error  within  pre¬ 
scribed  bounds.  This  method  requires  little  knowledge  of  the  problem 
being  solved  and  therefore  is  more  general  and  flexible.  This  last  approach 
was  the  one  pursued  in  the  study  to  be  described. 
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TRAJECTORY  PROBLEM 


The  trajectory  problem  can  be  divided  into  three  broad  cases. 

A  Six  degrees  of  freedom 
B  Five  degrees  of  freedom 
C  Three  degrees  of  freedom 

The  six  degrees  of  freedom  problem  in  its  most  general  form  is  a 
system  of  12  first  order  coupled  nonlinear  differential  equations.  Six 
of  these  equations  (3  second  order  D.  E. )  govern  the  motion  of  the 
center  of  gravity  and  the  remaining  six  govern  the  attitude  of  the  body. 
The  solutions  of  the  equations  governing  the  attitude  and  angular  rates 
are  oscillatory  in  nature,  the  oscillations  usually  having  small  periods 
in  the  independent  variable,  time.  These  small  periods  require  small 
mesh  size  to  trace  the  oscillation  and  to  maintain  stability.  This 
requirement,  coupled  with  the  usually  large  number  of  total  steps 
needed  to  take  the  trajectory  to  the  terminal  point  creates  significantly 
long  computation  times. 

The  major  difference  between  six  degree  of  freedom  trajectory 
problems  and  five  degree*  problems  is  the  choice  of  coordinate  systems 
specifying  the  orientation  of  the  body  axes.  For  spinning  missiles 
or  shell  tne  computational  problems  can  become  extremely  difficult 
when  the  kinematics  are  posed  in  the  six  degree  scheme.  The  attitude 
parameters,  for  example,  oscillate  sinusoidally  with  a  frequency  equal 
to  the  spin  rate.  This  spin  rate  can  often  be  much  higher  than  the  pitch 
and  yaw  frequencies  and  where  practical  the  problem  should  be  restated 
and  solved  using  five  degree  kinematics.  Both  of  these  simulations 
usually  include  additional  differential  equations  (in  addition  to  a  basic 
12  or  11)  for  representing  mass  changes,  guidance  and  control  etc. 

*• 

♦Five  degree  systems  have  also  been  called  '’fixed  plane  systems”8 
and  as  the  "Frick  slip  frame”9. 
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These  problems  are  usually  characterized  by  very  complicated  deriva¬ 
tive  sequences,  i.  e. ,  the  time  required  to  evaluate  the  derivative  is 
usually  large  in  comparison  to  the  time  required  to  execute  the  oper¬ 
ations  pertaining  to  the  integration  itself. 

The  three  degree  of  freedom*  problem  is  a  system  of  six  first  order 
ordinary  differential  equations  and  represents  an  approximation  to  the 
representations  described  above.  In  this  simulation  only  the  motion  of 
the  center  of  gravity  is  represented,  and  since  the  oscillations  are  not 
present  this  problem  is  not  nearly  so  demanding  of  computation  time. 

The  fact  that  very  large  numbers  of  these  trajectories  are  computed, 
however,  creates  a  need  for  economizing  the  solution  as  much  a n 
possible. 

METHODS  STUDIED 

Altogether,  eight  methods  were  studied  as  applied  to  trajectories, 
and  an  additional  method  (after  Scraton)  was  studied  on  a  prototype 
system  of  equations.  Six  of  these  methods  applied  to  trajectories  were 
predictor-corrector  type  schemes  and  among  these  six  were  some  which 
differed  only  in  the  sequencing  of  the  successive  application  of  the  pre¬ 
dictor  and  the  corrector  step;  the  other  two  methods  were  the  RKG 
method  and  the  Kutta-Merson  procedure. 


♦Recent  work  done  by  the  Firing  Tables  Branch  of  the  Computing 
Laboratory  has  been  successful  in  superimposing  upon  the  three  degree 
of  freedom  trajectory  for  spinning  shell  an  estimate  of  the  Myaw  of  repose' 
This  simulation  approaches  the  fidelity  of  the  six  and  five  degree  simu¬ 
lation  at  a  cost  not  much  greater  than  the  usual  three  degree  of  freedom 
simulation.  This  model  is  described  by  Lieske  and  Reiter10  . 


14 


DESCRIPTION  OF  THE  ALGORITHMS 


The  differential  equation  is  assumed  to  be  the  form 

Y'  =  f  (x,  Y) 


and  the  extension  to  a  system  is  done  in  the  usual  manner.  The  general 
formulas  for  the  predictor-corrector  scheme  are  as  follows: 

(P) 


=  ax 


r„_,  +  Y  +  dx  Y  +  h  (ex  f  +  fx  f 
n  1  n  c  n  j  n  n  1 


+  f„-2  +  kl  f„-3> 


(1) 


(c) 

Y„+l  ■  *•  Yn  +  b»  Yn-1  +  Yn“2  +  h  (e=  fn+l  +  f=  f„  +  8=  Vi 

+  ka  f  )  (2) 

n-c 

The  superscripts  (p)  and  (c)  denote  predicted  and  corrected  values 
respectively.  After  Hull  and  Creemer  we  adopt  the  notation  P  (EC)  m 
to  denote  the  mode  of  application  of  a  pair  of  formulas.  P  represents 
the  prediction  step,  E  the  evaluation  of  derivatives  and  C  the  corrector 
step,  m  is  the  number  of  times  the  cycle  (EC)  is  performed.  In  this 
study  two  values  of  m  were  utilized,  (m  =  1,  2).  The  RKGMA  procedure 
as  described  earlier  is  represented  by  the  cycle  PE  but  was  not  studied 
in  this  investigation  because  of  the  previous  observations  cited  and  it 
provides  no  estimate  of  the  truncation  error.  The  coefficients  for  the 
Adam's  formulas  in  equations  (1)  and  (2)  can  be  found  in  References  2 
and  3  and  the  Crane-Klopfenstein  coefficients  are  found  in  Reference  5. 
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The  Kutta-Merson  formulas  are  listed  below: 

KUTTA-MERSON  FORMULAS 
Yi  =  Y0  +  1/3  hf(x0,  Y0) 

Ya  =  Y0  +  1/6  h  f  (x0  ,  Y0  )  +  1/6  h  f  (x0  +  l/3h,  Yj. ) 

Y3  =  Yo  +  1/8  hf(x0,  Y0  )  +  3/8  h  f  (x0  +  l/3h,  Y2  ) 

Y4  =  Y0  +  1/2  hf(x0,  Y0  )  -  3/2  h  f  (x0  +  l/3h,  Ya  )  +2hf  (x0  +  l/2h,  Y3) 

Y6  =  Y0  +  l/6  h  f  (x0 ,  Y0  )  +  2/3  h  f  (x0  +  l/2h,  Y3)  +  l/6hf(x0+h,  Y4) 
The  procedure  provides  two  estimates  of  the  solution  at  the  end 
of  the  step,  namely  Y4  and  Yg.  If  the  differential  equation  is  of  the 


f°rm  Y*  =  f  (x,  Y)  = 

the  truncation  error  in  Y4  is 


a  x  +  b  Y  +  c 

i/l20  hB  Y^  and  in  YB  is 


1/720  h6  Y (gK 


The  use  of  this  error  estimate  for  systems  of  equations,  none  of  which 


have  the  above  form,  is  merely  a  computational  expediency  which  seems 
to  work  in  most  instances  and  should  accordingly  be  used  with  caution. 

LISTING  OF  METHODS 


Code 

Method 

Formulas 

1 

RKG 

Rung  e  -K  utta-G  ill 

2 

PEC 

Adams-Bashforth  Predictor 

3 

PECE 

Adams-Bashforth  Predictor 
Adams -Moulton  Corrector 

4 

PECEC 

Adams-Bashforth  Predictor 
Adams-Moulton  Corrector 

5 

PEC 

Crane-Klopfenstein  Predictor 

6 

PECE 

Crane-Klopfenstein  Predictor 
Adams-Moulton  Corrector 

7 

PECEC 

Crane-Klopfenstein  Predictor 
Adams-Moulton  Corrector 

8 

KM 

Kutta-Merson 
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ADJUSTMENT  OF  STEP  SIZE 

The  utilization  of  truncation  error  to  adjust  step  size  with  the 
predictor-corrector  methods  is  based  on  the  considerations  that  follow. 
See  Reference  2.  If  the  true  solution  at  the  point  n+1  is  denoted  as 

and  all  the  solution  and  derivative  values  occurring  in  the  right 
hand  members  of  (3)  and  (4)  were  known  exactly,  then 


(P) 


Y„+i  -  Y„+1  +  ^ 

|i  (c)  {P> 

Yn+1  =  Yn+ 1  +  E*  h6  +kh(fn+1-fn+i  ) 


(3) 

(4) 


k  is  equal  to  3/8  in  the  Adams  method  and  £  l  and  £a  both  lie  between 

x  „  and  x  ,  .  The  usual  procedure  is  to  assume  that  the  term 

n-3  n+1  * 

(c)  (p) 

Q  =  k  h  (f  .  -f  ,  ,  )  is  equal  to  zero,  a  condition  which  holds  only 
n+1  n+1  ^ 

after  iteration.  Since  the  iteration  requires  going  through  the  derivative 
evaluations  each  cycle  it  is  generally  more  efficient  to  proceed  with 
smaller  h  and  no  iteration.  The  term  Q  is  assumed  to  be  negligible. 

To  get  a  computable*  estimate  of  the  truncation  error  after  the  comple¬ 
tion  of  the  corrector  step  the  usual  assumption  is  that  h  is  sufficiently 
small  to  equate  Y^(£x  )  and  Y^  (£a  ) .  This  leads  to  the  equation 


hs  )  _  1 

5!  "  (Ei  -  Ea 


(c)  (p) 

(Yn+l  -  Yn+I> 


The  error  in  the  corrector  step  is  then  approximated  by 


T  « 


(Ex  -  Eg) 


(Y(‘> 

n+1  n+1 


(5) 


(6) 


♦The  ability  to  explicitly  evaluate  Y^  (£a )  is  not  practical  in  most 
cases  and  hence  is  not  considered  to  be  computable. 
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The  Kuita-Merson  procedure  is  not  really  a  predictor-corrector 
method  but  the  sequencing  of  operations  is  quite  similar,  i.  e. ,  the  two 
estimates  of  the  solution  at  n+1  are  obtained  and  the  difference  is 
utilized  in  estimating  the  truncation  error.  The  following  table  lists 
the  quantity  E2/(E1-E2)  for  the  various  methods. 


Method 

Adams 

C  rane-Klopfenstein 
Kutta-Merson 


E2/(E1-E2) 
1/4 
l/ 1 6 
1/5 


Truncation  Error 
19/720  h6  Y^(|) 
19/720  h5  Y®(i  ) 
1/720  hB 


The  classical  approach  to  step  size  adjustment  is  to  halve  or 
double  as  the  situation  warrants.  For  example,  if  T  denotes  the 
estimated  truncation  error  then  the  most  common  techniques  attempt 
to  bound  T  in  the  interval 


£  T  £  C 

h 


d  —  double 
h  —  halve 


(7) 


where  C,  denotes  the  critical  value  at  which  to  halve  and  C  ,  the  value 
h  d 

at  which  to  double.  The  motivation  for  doubling  and  halving  arises  from 
the  problems  caused  by  restarting.  For  instance  the  problem  of  re¬ 
starting  when  doubling  the  interval,  can  be  solved  by  storing  additional 
past  values  and  by  selecting  the  appropriate  values  for  the  new  step 
size  2h.  When  halving  one  can  use  interpolation  formulas  to  approximate 
the  required  data  for  the  new  interval  h/2.  Such  formulas  can  be  found 
in  Reference  11.  The  routines  utilized  in  this  study  were  started  by  the 
RKG  procedure  and  therefore  this  problem  of  restarting  was  not  present. 
Accordingly,  the  step  size  was  not  restricted  to  halving  and  doubling  but 
was  chosen  so  as  to  achieve  a  fixed  truncation  error.  The  step  size 
in  the  predictor-corrector  scheme  cannot  be  varied  continuously  because 
such  an  application  would  mean  that  most  of  the  computing  time  would  be 
spent  in  the  starting  phase.  For  the  predictor-corrector  scheme  the 
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step  size  was  changed  only  when  T  exceeded  the  bounds  of  (8) 


C  .  s;  T  ^  C, 
i  d 


1  —  increase 


d  —  decrease 


When  T  >  the  step  was  discarded  and  a  decreased  step  size  was 
computed.  When  T  £  C.  for  five  steps,  then  h  was  increased.  The  new  size 
was  always  chosen  so  as  to  move  T  back  to  the  middle  of  the  interval,  (8). 
The  KM  process  is  self  starting  and  hence  the  step  size  was  adjusted 
continuously.  The  new  step  size  is  determined  from  the  following  consid¬ 
eration.  The  quantity  Y^(£  )  is  assumed  to  vary  slowly  and  is  assumed 
constant  from  one  step  to  the  next.  The  old  step  size  is  denoted  as 

h0  and  the  new  as  h  .  T  is  the  maximum  truncation  error  occurring 

n  m 

for  those  equations  in  the  system,  h^  is  selected  to  make  the  truncation 

error  equal  to  T  . ,  the  desired  truncation  error.  T  and  T  ,  are  repre- 
^  d  m  d 

sented  by  1 


T  —  Eg  h  q 
m 

T,  =  Ea  h 


y@te ) 


Solving  (9)  and  (10)  for  h 


h  =  h0(T  /T  ) 
n  d  m 


Since  this  is  only  an  estimate  of  the  truncation  error  it  is  expected  that  at 

the  end  of  a  step  T  will  exceed  T,  as  often  as  it  is  less  than  T ,.  A 
r  m  d  d 

criterion  has  to  be  established  for  rejecting  and  repeating  a  step  when 

T  >  T,.  To  avoid  excessive  repeating  of  steps,  (11)  is  modified  as 
m  u 

follows:  . 

h  =  ho  (fT  /T  )  1/5  (12) 

n  am 

f  should  be  smaller  than  1  and  in  this  study  a  value  of  0.  1  was  used. 
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SETTING  THE  TRUNCATION  ERROR 
In  the  application  to  trajectories  the  control  on  the  step  size 
should  take  into  consideration  the  following  factors:  (1)  the  desired  end 
result  is  a  solution  giving  the  position  of  the  center  of  gravity  with  an 
error  no  greater  than  several  meters,  (2)  the  accuracy*  of  the  para¬ 
meters  governing  attitude  is  not  critical  since  the  attitude  merely  pro¬ 
duces  a  perturbation  or  second  order  effect  on  the  motion  of  the  center 
of  gravity.  Hence  the  consideration  here  is  to  set  the  truncation  error 
at  a  value  which  will  insure  stability;  i.  e. ,  high  accuracy  is  not  too 
important.  For  this  reason  the  algorithms  which  were  coded  for  the 
various  methods  were  designed  so  as  to  have  individual  controls  on  each 
equation.  This  was  executed  by  having  the  capability  of  specifying  a  T^ 
for  each  equation.  In  the  six  and  five-degree  study  the  controls  were 
effective  only  on  P,  Q,  R,  the  angular  rates  and  0,  6,  <f>  the  Euler  angles. 
See  appendix. 

DISCUSSION  OF  RESULTS 
Six  and  Five  Degrees  of  Freedom 

The  results  shown  in  Table  1  are  representative  of  the  relative 
economy  and  accuracy  of  the  eight  integration  methods  for  the  six  and 
five  degree  problems.  Trajectory  1  is  a  missile  trajectory  containing 
a  six  degree  phase  and  a  five  degree  phase.  The  running  time  spent 
in  the  boost  phase  (six  degrees)  is  approximately  the  difference  in  the 
running  time  between  trajectory  I  and  trajectory  II.  Trajectory  II 
contains  only  a  five  degree  phase.  Integration  type  1  was  the  RKG  method 
with  fixed  step  size  and  was  used  as  the  standard  of  comparison.  For 


*In  the  absence  of  a  true  solution  the  accuracy  can  be  defined  by  comparing 
results  against  a  standard  solution.  If  a  given  problem  is  integrated  using 
a  particular  method  for  several  step  sizes,  then  for  a  region  of  h  the 
solution  will  not  vary  appreciably.  In  this  region  one  can  choose  the 
standard  value. 
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integration  type  1  the  integration  interval  was  held  constant  at  0.  01 

AT 

during  the  six  degree  boost  phase  and  at  0. 1  during  the  five  degree  phase 
thereafter.  The  results  could  be  summarized  as  follows.  Of  the  pre¬ 
dictor-corrector  methods  listed  the  results  are  not  very  conclusive 
except  to  indicate  that  a  method  using  two  evaluations  is  probably  better 
than  a  method  using  only  one  evaluation.  This  conclusion  is  drawn  from 
the  observation  that  the  running  time  and  the  accuracy  are  not  signifi¬ 
cantly  different  and  two  evaluations  would  generally  be  better,  due  to 
considerations  of  safety  arising  from  greater  stability.  No  further 
attempt  was  made  to  choose  between,  the  various  predictor-corrector 
methods  because  of  the  fact  that  KM  by  far  seemed  the  best.  This  is 

especially  true  when  the  accuracy  requirement  is  made  less  stringent. 

-3 

For  truncation  errors  of  10  ,  KM  is  almost  twice  as  efficient  as  any 

predictor-corrector  method. 

Figure  1  shows  the  step  size  history  along  trajectory  II,  using  the 

KM  procedure  and  the  Adams  predictor-corrector  method,  PECE.  In 

-3 

both  instances  the  truncation  error  was  set  to  10  .  The  data  for  KM 

is  plotted  only  at  intervals  of  five  seconds. 
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TABLE  1 


TRAJECTORY  I 

SIX  ANO 

FIVE  DEGREES 

INTEGRATION 

TYPE 

TRUNCATION 

ERROR 

TIME  OF 
FLIGHT  ERROR 
SECONDS 

DEFLECTION 

ERROR 

METERS 

RANGE 

ERROR 

METERS 

COMPUTATION 

TIME 

MIN 

1 

2 

3 

4 

5 

6 

7 

8 

10  ~4 

I  t 

•  • 

•  1 

•  t 

•  • 

•  • 

.0001 

-.0009. 

.0008 

.0001 

-.0010 

.0009 

.0001 

.1 

1.0 

.4 

.0 

1.1 

.9 

.0 

.1 

.2 

.1 

.1 

.2 

-.1 

.0 

4.16 

2.03 

1.85 

2.01 

2.08 

2.04 

2.01 

1.76 

I 

m 

4.16 

2 

10  “3 

.0034 

2.9 

2.1 

1.77 

3 

•  • 

-.00A9 

3.7 

1.8 

1.49 

A 

•  • 

.0032 

.4 

.2 

1.71 

5 

•  t 

.0014 

2.9 

2.8 

1.83 

6 

«  • 

-.0044 

4.0 

1.7 

1.72 

7 

1  • 

.0014 

1.1 

.6 

1.64 

8 

t  f 

-.0013 

.5 

_ 

.2 

.93 

SIX  DEGREES 


-4 

5.78 

8 

10  , 

-.0003 

3.4 

.0 

8 

10 

-.0249 

22.5 

-2.2 

3.28 

THE  TRAJECTORY  MODEL  ON  WHICH  THIS  DATA  WAS  COLLECTED  WAS  THAT 
OF  THE  LANCE  MISSILE  SYSTEM.  TRAJECTORY  l,  WHEN  USING  THE  RKG 
INTEGRATION  ANO  SIX  DEGREES  OF  FREEDOM  WITH  A  CONSTANT  STEP  SIZE 
OF  .025  SECONDS  REQUIRED  APPROXIMATELY  15  MINUTES  SOLUTION  TIME. 
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TABLE  I 
(CONT.  J 


INTEGRATION 

TRUNCATION 

TRAJECTORY  II 

FIVE  DEGREES 

TIME  OF  DEFLECTION 

RANGE 

COMPUTATION 

TYPE 

ERROR 

FLIGHT  ERROR 

ERROR 

ERROR 

TIME 

SECONDS 

METERS 

METERS 

MIN 

1 

i"  " 

~  -4 

! 

2.84 

2 

10 

.2 

.1 

1.57 

|  3 

•  • 

.0009 

1.1 

.2 

1.34 

i  4 

1  « 

.0008 

.5 

.1 

1,50 

!  5 

1  • 

-.0001 

•1 

.1 

1.59 

6 

•  1 

-.0010 

1.2 

.2 

1.42 

7 

.0009 

.9 

-.1 

1.40 

8 

•  I 

.1 

.0 

1.13 

1 

!  I 

■  -3 

2.84 

:  2 

10 

.0041 

3.3 

2.1 

1.4  6 

:  3 

f  • 

-.0045 

3.9 

1.8 

1.07 

.0031 

.4 

.2 

1.27 

5 

f  • 

.0022 

3.0 

3.5 

1.53 

j  6 

•  1 

-.0043 

4.0 

1.33 

1  7 

•  1 

.0014 

.2 

1.33 

1  8 

•  • 

.0016 

.5 

.72 

SIX  DEGREES 


H 

- =3 - 

10-3 

3.6 

5.17 

j 

8 

1 _ 

10 

■Bin 

21.8 

MEM 

2.96 
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STEPSIZE  h 
(SECONDS) 


5°  TRAJECTORY 


.  KUTTA  MERSON 
-  METHOD  3 


10  ?0  30  40  50  60  70  80  90  100  110  120 

TIME  ALONG  TRAJECTORY 


(SECONDS) 
FIGURE  1 
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Three  Degrees 

As  pointed  out  previously,  the  most  critical  running  time  problem 
occurs  in  the  five  and  six  degree  trajectory  models.  Accordingly  the 
emphasis  was  placed  on  these  problems  initially  and  after  KM  was 
selected  as  the  most  efficient  method,  it  was  applied  to  the  three  degree 
problem.  The  hope  was  that  Kutta-Merson  would  prove  to  be  at  least  as 
efficient  as  previously  used  methods  for  three  degrees  thus  allowing  the 
use  of  a  single  method  for  all  trajectory  programs.  The  equations 
for  the  model  used  in  this  study  are  found  in  Reference  10.  In  this  study 
KM  was  compared  against  a  method  currently  used  in  the  Firing  Tables 
Branch  for  the  mass  production  of  all  firing  tables  for  cannon  artillery. 
This  method  is  not  widely  used  elsewhere  ,  but  was  tailored  to  this  parti¬ 
cular  application.  The  method  is  a  self  starting  predictor-corrector 
method.  The  three  degree  equations  are  of  the  form 

•  •  ^  ^  m  • 

X  =  F  (  X,  X  )  where 

v 

X  is  the  acceleration  vector ,  X  the  velocity  and  X  is  the  position  vector. 
The  integration  algorithm  is  as  follows. 


— ♦ 

X 

— * 

xn 

-4 

+  h  Xft 

p 

• 

0 

—4 

X 

V 

+  h  Xrt 

p 
•  • 

0 

0 

.  • 

-4 

X 

= 

F  ( 

X  ,  X  ) 

p 

p  p 

•  • 

-♦ 

X 

c 

• 

-4 

X 

p 

• 

-h/2  (X0  -  Xp  )  + 

#  •  •  •  # 

-4 

X 

c 

•  • 

= 

-4 

X 

p 

-h/2  (  XQ  -  5p  ) 

• 

—4 

X 

c 

= 

F( 

X  ,  X  ) 
c  c 
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The  quantities  with  the  subscript  0  denote  values  at  the  beginning  of  a  step. 
The  subscript  p  denotes  predicted  values  and  c  denotes  corrected  values. 

The  predictor  formula  is  that  of  Euler  and  has  a  truncation  error  of  order 

2  -*  5 

h  .  The  corrector  formula  for  X  is  of  order  h  and  that  of  X  of  order 

2  C  C 
h  .  The  fact  that  the  predictor  and  corrector  formulas  for  X  are  of 

different  order  precludes  the  estimation  of  truncation  error  by  the  usual 

methods.  Previous  techniques  for  adjusting  step  size  with  this  method 

were  based  on  an  empirical  formula  which  expressed  h  as  a  function  of 

Mach  number.  Each  weapon  had  its  own  function,  a  typical  one  being  as 

follows: 


h 

s 

.  5 

0.0 

£ 

M 

< 

.764 

h 

s 

.25 

0.  764 

£ 

M 

< 

.936 

h 

s 

.25 

0.936 

£ 

M 

< 

.960 

h 

= 

.0625 

0.960 

£ 

M 

< 

.984 

h 

s 

.0625 

0.984 

£ 

M 

< 

1.016 

h 

s 

.  125 

1.016 

£ 

M 

< 

1.  25 

h 

- 

.  25 

1.25 

£ 

M 

< 

4.0 

M  =  Mach  Number 

The  function  above  was  used  in  the  three  degree  study  for  the  175mm  Gun 
and  a  similar  function  was  used  for  the  155mm  Gun.  The  integration 
technique  described  above  with  step  size  adjusted  as  a  function  of  Mach 
number  will  be  referred  to  as  method  A. 

An  alternative  procedure,  which  has  been  used  on  some  problems 

^4  “4 

in  the  Computing  Laboratory,  is  to  use  the  formula  h  =C  I  X  /  X  | 
to  continually  adjust  h.  The  constant  C  is  empirical  and  is  determined 
by  experiment.  This  formula  and  the  integration  technique  described 
above  will  be  referred  to  as  method  B. 

The  three  degree  study  consisted  of  three  methods.  The  first  method 
was  the  Kutta-Merson  procedure,  the  second  was  method  A  and  the  third 
was  method  B.  Initially,  experiments  were  conducted  to  determine  control 
values  for  the  truncation  error  which  would  consistently  yield  good 


accuracy  without  sacrificing  excessive  computation  time.  The  result  of 
these  experiments  indicated  that  a  good  choice  of  truncation  error  with 
the  Kutta  Merson  procedure  was  0.002  on  the  three  differential  equations 
associated  with  the  velocity  and  0.  1  on  the  three  equations  associated 
with  the  position.  For  Method  B  a  value  of  C  =  0.02  was  selected.  The 
three  methods  were  compared  on  two  weapons,  the  175mm  Gun  and  the 
1 55mm  Gun.  Three  velocity  levels  and  three  quadrant  elevations  were 
studied  for  a  total  of  18  trajectories  for  each  of  the  three  methods.  Each 
trajectory  was  run  50  times  so  that  a  reliable  estimate  of  the  computation 
time  would  be  obtained.  The  results  are  listed  in  Tables  2  and  3  and  can 
be  summarized  as  follows: 

For  the  155mm  Gun,  Method  A  seems  to  be  as  fast  as  Kutta-Merson 
or  Method  B  but  the  accuracy  is  not  acceptable.  To  reduce  the  error 
using  Method  A  would  require  longer  solution  time.  This  is  evidenced 
by  the  data  in  Table  3  for  the  175mm  Gun  where  the  accuracy  for  Method  A 
is  acceptable  but  the  solution  time  is  not  competitive  with  respect  to  the 
Kutta-Merson  procedure  or  Method  B.  In  comparing  Method  B  against 
Kutta-Merson  it  seems  that  no  significant  difference  in  accuracy  or  speed 
is  apparent. 

The  results  of  this  study  indicate  that  Kutta-Merson  is  a  better 
procedure  for  the  three  degree  computations  than  the  technique  presently 
used;  i.  e. ,  the  adjustment  of  h  as  a  function  of  Mach  number,  and  at 
least  as  efficient  and  accurate  as  Method  B  described  earlier.  Since 
it  is  desirable  to  standardize  on  a  method  for  various  problems  within 
an  organization,  Kutta-Merson  by  virtue  of  being  more  efficient  for  six 
degrees  and  at  least  as  efficient  for  three  degrees  would  seemingly  be  the 
best  method.  Another  important  point  at  least  in  the  BRL  is  the  following: 
In  some  rocket  trajectories  the  simulation  is  done  in  stages,  five  degrees 
for  the  first  phase,  (boost)  and  three  degrees  for  the  second  phase  (coast). 
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Kutta-Merson  can  handle  both  stages  with  equal  ease. 

Figure  2  is  a  trace  of  the  step  size  variation  as  a  function  of  the 
time  along  the  trajectory  for  the  155mm  Gun.  This  trajectory  had  an 
elevation  of  1124  mils  and  a  muzzle  velocity  of  561  meters/sec.  The 
truncation  error  was  set  at  0.01  on  all  six  differential  equations  of  the 
three  degree  system. 
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TABLE  2 


THREE  OEGREES 


WEAPON 

MUZZLE 

VEL 

M/S 

QUAD 

ELEV 

MILS 

COMP  TIME 
FOR 

50  TRAJ 

MIN 

TIME  OF  TIME  OF 
FLIGHT  FLIGHT 
ERROR 

SECONDS  SECONDS 

KUTTA-MERSON 

RANGE 

M 

RANGE 

ERROR 

M 

DEFL 

M 

DEFL 

ERROR 

M 

155 

280.0 

200 

.38 

11.015 

.000 

2870.05 

.01 

11.24 

.04 

155 

280.0 

800 

.59 

38.557 

.001 

6684.86 

-  .09 

144.31 

.67 

155 

280.0 

1200 

1.05 

50.116 

-  .001 

4541.49 

-  .95 

270.92 

.87 

155 

378.0 

200 

.69 

13.965 

.000 

4435.45 

.03 

21.40 

.06 

155 

378.0 

300 

.99 

48.002 

.001 

10016.63 

-  .07 

237.25 

1.02 

155 

378.0 

1200 

1.52 

62.594 

.000 

6916.38 

-  1.03 

441.10 

1.25 

155 

684.3 

200 

.79 

22.797 

.001 

9764.53 

.29 

75.98 

.91 

155 

684.3 

800 

1.55 

68.100 

.001 

18020.35 

-  .11 

635.09 

3.04 

155 

684.3 

1200 

2.73 

90.005 

.000 

12530.35 

-  3.07 

1093.87 

2.74 

T0TAL=10.29 


METHOD  A 


155 

280.0 

— 

.22 

11.010 

— 

.004 

2868.74 

-1.24 

11.01 

-  .10 

155 

280.0 

.49 

38.551 

- 

.005 

6683.93 

-  1.01 

143.72 

.26 

155 

280.0 

.60 

50.112 

- 

.005 

4539.97 

-  2.37 

270.96 

.98 

155 

378.0 

.83 

13.964 

- 

.001 

4435.17 

-  .15 

21.05 

-  .17 

155 

378.0 

.91 

47.998 

- 

.003 

10016.33 

-  .42 

236.41 

.39 

155 

378.0 

1.05 

62.592 

- 

.002 

6915.09 

-  2.26 

441.18 

1.40 

155 

684.3 

L.00 

22.793 

- 

.002 

9762.55 

-  1.57 

75.36 

.37 

155 

684.3 

800 

3.11 

68.097 

- 

.007 

18019.66 

-  .76 

633.95 

1.88 

155 

684.3 

1200 

2.94 

90.005 

.000 

12520.28 

-12.35 

1095.25 

3.95 

T0T4L= 11 .15 


METHOD  6 


155 

280.0 

200 

.28 

11.012 

— 

.002 

— 

2869.32 

"" 

.66 

10.53 

-  .58 

155 

280.0 

800 

1 

.00 

38.555 

- 

.001 

6684.79 

- 

.15 

142.56 

-  .90 

155 

280.0 

1200 

1 

.78 

50.116 

- 

.001 

4542.19 

- 

.15 

269.42 

-  .56 

155 

378.0 

200 

.32 

13.961 

- 

.004 

4433.93 

- 

1.39 

20.73 

-  .49 

155 

378.0 

800 

1 

.09 

48.000 

- 

.001 

10016.77 

.02 

235.49 

.53 

155 

378.0 

1200 

1 

.84 

62.592 

- 

.002 

6916.90 

- 

.45 

439.60 

-  .18 

155 

684.3 

200 

.58 

22.792 

- 

.003 

9762.33 

- 

1.79 

74.85 

-  .14 

155 

684.3 

800 

1 

.34 

68.095 

- 

.003 

18019.48 

- 

.94 

633.19 

1.12 

155 

684.3 

1200 

2 

.05 

89.998 

— 

.007 

12529.27 

3.36 

1092.60 

L 

1.30 

T0TAL=10.28 
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TABLE  3 


THREE  OEGREES 


WEAPON  MUZZLE 
VEL 

QUAD 

ELEV 

COMP  TIME 
FOR 

50  TRAJ 

TIME  OF 
FLIGHT 

TIME  OF 
FLI'HT 
ERRv.i 

RANGE 

RANGE 

ERROR 

OEFL 

DEFL 

ERROR 

H/S 

MILS 

MIN 

SECONDS 

SECONDS 

M 

M 

M 

H 

KUTTA-MERSON 


175 

510.5 

200 

.70 

18.643 

7248.86 

.35 

25.76 

.43 

175 

510.5 

800 

1.27 

59.722 

15186.80 

.05 

194.91 

.57 

175 

510.5 

1200 

1.60 

75.540 

1  1* 

10999.15 

.51 

388.68 

-  .85 

175 

704.1 

20-., 

.74 

24.788 

if] 

11852.43 

.88 

64.98 

1.26 

175 

70  4.1 

300 

1.53 

76.577 

22212.09 

-  .05 

410.69 

3.34 

175 

704.1 

1200 

2.00 

100.920 

.001 

16738.79 

.86 

697.62 

.99 

175 

914.4 

200 

.82 

30.897 

.004 

17243.61 

1.43 

122.26 

2.25 

175 

914.4 

800 

1.66 

96.351 

.007 

32223.58 

.31 

882.43 

-1.50 

175 

914.4 

1200 

2.13 

127.643 

.000 

26104.62 

1.99 

1291.52 

3.86 

T0TAL=12.45 


METHOD  A 


175 

mm 

200 

.90 

18.640 

,  _ 

.002 

7247.68 

.65 

25.04 

-  .12 

175 

■SIS 

800 

1.81 

59.718 

- 

.003 

15185.48 

- 

1.24 

193.56 

-  .12 

175 

510.5 

1200 

1.98 

78.534 

- 

.007 

10997.56 

- 

1.16 

387.98 

-1.22 

175 

704.1 

200 

.81 

24.783 

- 

.002 

11850.24 

- 

1.12 

63.77 

.24 

175 

n 

800 

3.38 

76.572 

- 

.003 

22210.82 

- 

1.28 

407.95 

.85 

175 

704.1 

1200 

2.80 

100.913 

- 

.006 

16737.49 

- 

.52 

695.63 

-  .70 

175 

914.4 

200 

.89 

30.892 

- 

.001 

17240.73 

- 

1.29 

120.68 

.84 

175 

914.4 

800 

3.43 

96.342 

- 

.003 

32220.39 

- 

2.77 

876.99 

3.02 

175 

914.4 

1200 

3.18 

127.632 

— 

.010 

26104.50 

1.63 

1287.88 

.46 

TOTAL=19.18 


METHOD  B 


175 

510.5 

mm 

.41 

18.636 

— 

.006 

7246.09 

2.24 

24.19 

-  .97 

175 

510.5 

1.16 

59.717 

- 

.004 

15186.23 

- 

.49 

192.26 

-1.71 

175 

510.5 

RES 

1.87 

78.533 

- 

.008 

10998.49 

- 

.23 

387.82 

-1.38 

175 

704.1 

200 

.54 

24.780 

- 

.005 

11848.74 

- 

2.62 

62.58 

-  .95 

175 

704.1 

800 

1.30 

76.569 

- 

.006 

22210.49 

- 

1.61 

406.82 

-  .28 

175 

704.1 

1200 

1.94 

100.904 

- 

.015 

16737.25 

- 

.76 

695.73 

-  .60 

175 

914.4 

200 

.65 

30.888 

- 

.005 

17239.26 

- 

2.76 

119.38 

-  .46 

175 

914.4 

800 

1.41 

96.335 

- 

.010 

32219.31 

- 

3.85 

877.87 

3.90 

175 

. 

914.4 

1200 

2.01 

127.617 

— 

.025 

26100.67 

— 

2.20 

1288.65 

1.23 

T0TAL=11.29 
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EFFICIENCY  OF  THE  METHODS 


Merson6  observed  that  for  dynamic  problems,  one  could  use  with 
the  Runge-Kutta  method,  an  integration  interval  perhaps  seven  times 
larger  than  with  the  Milne  predictor.  This  general  observation  seems  to 
be  borne  out  by  this  study. 

The  results  are  somewhat  surprising  in  that  it  is  common  to 
attribute  to  the  predictor-corrector  methods  a  greater  efficiency  than 
the  Runge-Kutta  methods.  One  usually  assumes  that  when  the  order  of 
two  methods  is  the  same  and  one  has  say  two  evaluations  in  one  method 
and  four  in  the  other  then  the  method  with  two  evaluations  on  the  surface 
is  twice  as  efficient  as  the  other.  Hence  one  is  inclined  to  attribute  the 
contrary  results  of  this  study  to  considerations  of  stability.  This  is 
probably  true  but  two  other  factors  seem  important. 

The  first  of  these  is  that  to  avoid  the  cost  of  restarting  with  the 
predictor-corrector  methods,  one  has  to  employ  the  bounding  conditions 
specified  in  equation  (8).  This  means  that  one  is  not  always  using  the 
largest  step  size  consistent  with  the  desired  truncation  error  as,  e.g. , 
one  uses  with  the  KM  procedure. 

The  second  factor  is  accuracy  itself.  If  RKGMA,  the  conventional 
Adams  method,  and  KM  produced  the  same  accuracy  for  a  given  step 
size  then  one  might  choose  to  assign  to  each  method  a  theoretical  cost 
index  as  shown  in  the  second  column  of  Table  3. 

TABLE  3 


Method 

RKGMA 

Adams 

Kutta-Merson 


Cost  Index 
1.0 
2.0 
5.0 


Cost  Accuracy  Index 
1.0 
1.  2 
1.  7 
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Since  the  methods  do  not  yield  the  same  accuracy  (despite  all  being 
fifth  order)  for  the  same  step  size,  it  is  more  realistic  to  define  a  cost- 
accuracy  index.  A  cost-accuracy  index  can  be  defined  as  the  quantity 
N/R  where  N  is  the  number  of  derivative  evaluations  per  step  and  R  is  a 
step  size  ratio  required  to  achieve  a  fixed  truncation  error.  This  ratio 
can  best  be  found  by  comparing  methods  against  a  standard  method.  For 
the  purpose  of  comparison  let  the  RKGMA  method  be  the  standard  with  a 
cost  accuracy  index  of  1.0.  KM  has  an  N  of  5  and  can  allow  a  step  size 
l/5 

of  (251)  «  3  times  larger  than  RKGMA  to  achieve  the  same  truncation 

error.  This  ratio  is  obtained  by  equating  the  truncation  error  of  the 
two  methods  and  solving  for  the  step  size  ratio.  The  cost-accuracy 
index  for  the  three  methods  is  shown  in  the  third  column  of  Table  3. 

Even  with  the  more  realistic  comparison  of  a  cost-accuracy  index 
the  Kutta-Merson  procedure  shows  up  as  the  least  efficient  method.  The 
contrary  results  of  this  study  show  that  for  dynamic  problems,  such  as 
those  studied  here,  efficiency  can  only  be  determined  by  experimentation. 

The  cost-accuracy  index,  however,  suggests  that  the  additional 

evaluations  do  more  than  contribute  to  stability,  they  also  increase 

accuracy,  and  therefore  efficiency.  An  additional  factor  which  should 

be  noted  is  the  following:  In  representing  the  propulsion  and  aerodynamic 

characteristics  of  a  missile  ,  one  is  usually  faced  with  the  problem  of 

approximating  the  data,  by  some  curve  fitting  technique  so  that  it  can  be 

utilized  in  the  trajectory  problem.  Because  of  the  severe  nonlinearities 

of  the  data  and  because  of  practical  considerations  it  is  customary  to 

approximate  the  data  by  piecewise  polynomials.  At  the  juncture  of  these 

polynomials  the  derivatives  are  not  continuous  and  hence  one  of  the  basic 
* 

assumptions  of  any  integration  technique,  i.  e. ,  the  continuity  of  higher 
derivatives,  is  violated.  Intuitively  it  would  seem  that  an  integration 
method  which  uses  information  only  within  one  interval  (Runge-Kutta) 
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as  compared  to  one  which  uses  information  over  several  intervals 
(predictor-corrector)  would  be  less  affected  by  this  problem.  In  this 
study,  the  data  in  the  six  and  five  degree  problems  were  represented  by 
piecewise  linear  segments.  In  the  three  degree  problem  piecewise  4th 
degree  polynomials  were  used. 

Since  the  completion  of  this  study  the  author  has  become  aware  of 
work  done  by  Scraton7  who  obtains  an  estimate  of  the  truncation  error 
for  the  Runge-Kutta  method  when  making  the  same  assumptions  as  does 
Merson. 

This  assumption  is  that  the  differential  equation  is  the  form 

Y1  =  ax  +  bY  +  c  (13) 

where  a,  b,  and  c  are  constants. 

In  addition  Scraton  has  derived  Runge-Kutta  type  formulas  which  require 
five  evaluations  and  which  have  a  computable  estimate  of  the  truncation 
error  not  based  on  the  above  assumptions,  i.  e. ,  they  have  general 
validity.  This  recent  work  raises  the  question  as  to  whether  the  Runge- 
Kutta  procedure  with  Scraton' s  truncation  error  estimate  might  not  be 
a  better  method  than  the  KM  procedure  since  it  requires  only  four 
derivative  evaluations.  Furthermore,  might  not  Scraton' s  new  formulas 
be  more  attractive  because  of  their  general  validity.  Mr.  Glenn  Beck  cf 
the  Computing  Laboratory  has  observed  that  the  KM  procedure  as  used 
with  a  prototype  differential  equation  yields  greater  accuracy  than  does 
RK  or  Scraton' s  method  and  the  difference  in  accuracy  compensates  for 
the  fifth  evaluation.  The  belief  is  e.  g. ,  see  Fox13 ,  that  when  the 
system  of  equations  does  not  follow  the  form  (13),  that  Merson' s 
truncation  error  estimate  is  an  overestimate,  i.  e. ,  it  is  safe. 
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CONCLUSIONS 


In  five  and  six  degree  of  freedom  problems  the  Kutta-Merson 
integration  procedure  as  compared  against  the  fifth  order  predictor- 
corrector  techniques  seems  to  be  twice  as  economical  in  terms  of 
computation  time  when  applied  in  a  manner  which  allows  continuous 
interval  adjustment.  When  compared  against  the  RKG  procedure  with 
constant  step  size,  the  data  for  this  report  indicate  savings  of  almost 
4  to  1.  In  most  cases  the  ratio  would  be  higher  because  one  usually 
cannot  justify  an  extensive  stability  study  to  find  the  largest  value  of 
h  required  to  yield  the  desired  accuracy  and  maintain  stability. 

For  three  degree  problems  the  Kutta-Merson  procedure  seems 
to  be  more  efficient  than  presently  used  techniques  and  at  least  as 
efficient  as  a  proposed  modification. 

The  adoption  of  Kutta-Merson  as  the  standard  subroutine  for  all 
trajectory  computations  can  result  in  significant  savings  of  computation 
time  in  addition  to  the  benefits  gained  from  the  standardization.  In 
addition  it  has  already  proved  quite  useful  on  a  number  of  other  problems. 
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This  study  has  resulted  in  the  Kutta-Merson  nrocedure  being  coded 
as  an  integration  package  which  is  now  part  of  the  BRLESC  compiler.  The 
subroutine  description  can  be  obtained  by  contacting  the  Computation 
Branch  of  the  Computing  Laboratory. 
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APPENDIX 

EQUATIONS  OF  MOTION 
List  of  Symbols 
Position  in  inertial  space 

•4-4-4 

Propulsion  forces  about  body  axes  x*  y,  z 

•4  —4  -4 

Aerodynamic  forces  about  body  axes  x,  y,  z 


Components  of  gravitational  vector  along 

— 4  *4  ^ 

inertial  axes  x0  »  yo  >  z0 

—4  —4  -4 

Propulsion  moments  about  body  axes  x,  y,  z 

-4-4-4 

Aerodynamic  moments  about  body  axes  x,  y,  z 


Angular  rates  of  missile  about  body  axes,  x,  y,  z 

—4  «— 4  —4 

Moments  of  inertia  about  body  axes  x,  y,  z 


Euler  angles 

Transformation  matrix  between  body  and  inertial 
axes 

Unit  vectors  defining  the  inertial  axes 
Unit  vectors  defining  the  body  axes 
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LISTING  OF  EQUATIONS 
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The  equations  of  motion  listed  above  are  those  used  in  the  six  and 
five  degree  study.  These  equations  are  general  and  can  be  used  for  any 
missile  or  projectile  which  has  axial  symmetry.  The  detailed  listing 
of  the  computations  involved  in  obtaining  the  propulsion,  aerodynamic 
and  gravitational  forces  and  moments  are  not  included,  but  as  pointed 
out  earlier  these  computations  are  very  lengthy.  The  transformation 
matrix  relating  the  body  axes  and  the  inertial  axes  is  obtained  by  an 
Euler  angle  system.  The  equations  are  posed  for  programming  conven¬ 
ience  so  that  the  program  can  switch  from  six  degrees  to  five  degrees 
with  only  minor  changes  in  program  logic. 

These  equations  are  conventional  except  for  the  presence  of  P 
in  the  last  two  equations  of  (A2).  For  the  six  degree  system  P  is 
equal  to  P .  For  the  five  degree  system  P  is  determined  by 

P  =  -  R  tan  0  (A5) 

In  addition,  for  the  five  degree  system  <f>  =  ’<f>  =  0  in  (A3)  and  (A4). 

The  singularity  which  can  occur  in  the  $  equation  when  0  approaches 

90°  is  avoided  by  utilizing  a  "switch  over"  set  of  intermediate  axes 

0  *  0  £  90°. 

c 

It  should  be  noted  that  alternative  techniques  for  obtaining  [C  ]  are 
often  used.  A  method  which  seems  to  be  used  as  often  as  the  Euler 
angle  method  is  to  integrate  the  direction  cosine  derivatives  directly  by 
the  equations 
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This  method  avoids  the  evaluation  of  the  trigonometric  terms  in.[C] 
at  the  expense  of  the  additional  differential  equations.  The  experience 
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of  various  Government  agencies  and  the  industry  with  this  technique 
however  is  that  [C  ]  tends  to  gradually  deviate  from  orthogonality 
(especially  for  spinning  missiles)  and  accentuates  instability.  This 
problem  is  avoided  by  inserting  into  the  program  an  orthogonality 
loop  such  as  that  described  in  Reference  15  which  at  each  step  of  the 
integration  forces  [C]  back  to  a  condition  of  orthogonality.  Equations 
(A6)  can  be  utilized  with  both  five  and  six  degree  systems.  For  five  degree 
systems  P  is  found  from  the  equation 

P  =  -  R  tan  1 3/  m3  (A7) 

An  interesting  observation  is  that  conventional  six  degree  systems  such  as 
that  described  by  Zaroodney11  which  utilize  (A6)  for  finding  [C  ]  can  be 
modified  to  five  degrees  by  simply  replacing  P  with  P  as  noted  in 
(A6)  and  defining  P  as  in  (A7). 

The  rigid  body  system  of  equations  most  widely  used  in  firing 

1  a 

table  computations  at  BRL  is  that  described  by  Lieske  and  McCoy 
and  due  originally  to  H.  L.  Reed.  This  system  is  a  five  degree  one 
and  integrates  the  direction  cosine  derivatives  but  does  not  employ  the 
"body  axes"  as  is  usually  done,  i.  e.,  all  integrations  are  done  in  inertial 
space.  The  method  enjoys  the  advantage  of  being  simple  but  suffers 
the  disadvantage  of  being  unconventional.  The  missile  engineer  with 
the  customary  education  in  the  aeronautical  schools  usually  prefers  to 
express  the  forces  and  moments  in  terms  of  parameters  associated 
with  the  "body  axes"  such  as  P,  Q,  R,  the  angular  rates,  a,  0,  the 
angle  of  attack  and  angle  of  sideslip  and  0  ,  6  ,  <p ,  the  Euler  angles. 

This  seems  especially  true  for  guided  missile  systems  which  have 
"on  board"  inertial  platforms  which  in  some  way  or  other  are  "tied" 
to  the  body  axes. 
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